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EXACT SIMULATION OF DIFFUSIONS 

By Alexandros Beskos^ and Gareth O. Roberts 

Lancaster University 

We describe a new, surprisingly simple algorithm, that simulates 
exact sample paths of a class of stochastic differential equations. It 
involves rejection sampling and, when applicable, returns the location 
of the path at a random collection of time instances. The path can 
then be completed without further reference to the dynamics of the 
target process. 

1. Introduction. Exact simulation of stochastic differential equations 
(SDEs) is a notorious problem within the applied probability community. 
The objective of this paper is to present a first step toward the solution of 
this problem. We describe a new algorithm we call the Exact Algorithm for 
simulating a class of SDEs. It involves rejection sampling and, when appli- 
cable, returns exact draws from any finite-dimensional distribution of the 
solution of the SDE. 

Let B = {Bt; <t <T} be a scalar Brownian motion. Consider the gen- 
eral type of the one-dimensional Ito diffusion: 

(1) dXt = h{Xt) dt + a{Xt) dBt, <t <T,Xo = x eR 

for drift coefficient 6 : R i— > R and diffusion coefficient o" : R ^ R. Under cer- 
tain regularity conditions on b and a it can be shown that (1) has a solution 
{Xt',0 <t<T} weakly unique, that is, all the solutions have identical finite- 
dimensional distributions. Weak uniqueness, relatively more general than 
pathwise uniqueness, is sufficient for simulation purposes. For a formal def- 
inition of (1) see, for instance, [7]. 

Mathematical models of this kind are used to describe the evolution of 
stochastic phenomena in a wide range of disciplines and most times it is 
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important to be able to simulate them. In rare cases (1) has an explicit so- 
lution with identifiable transition density. However, in most of the models 
used in practice it is necessary to resort to a numerical solution of (1). This 
has traditionally implied the use of some of the time discrete approximation 
methods (Euler, Taylor's expansion, etc.) which rely on small time approx- 
imate increment distributions for the diffusion (for a detailed account of 
these methods see [5]). In many cases throughout this paper the results of 
our algorithm will be compared with the Euler scheme which approximates 
(1) via the recursion 



where we denote by AA(/i, S) the Gaussian distribution with mean /i and 
variance S. 

For a large class of processes the Exact Algorithm provides an alternative 
which involves no approximation (apart from that inherent in any computer 
simulation) and yet is computationally highly efficient. It returns skeletons 
of exact paths which can be easily filled in without further reference to the 
diffusion dynamics. Thus the method can be used to simulate the diffusion 
at a prescribed collection of time points, or alternatively at times which 
occur to be interesting to the user, after the completion of the algorithm. 

We begin (Section 2) by stating the rejection sampling technique in a 
way that serves our purposes. In Section 3 we present our method and in 
Section 4 give some results related with its efficiency. In Section 5 we apply 
the algorithm to a specific SDE and in Section 6 we take advantage of the 
properties of the Exact Algorithm to simulate exactly extremes and hitting 
times for the same SDE. Finally (Section 7), we present some ideas that 
could, in future research, overcome the restrictions of the algorithm and 
give some general conclusions. We are going to restrict the exposition of our 
algorithm to the case when cr = 1 . This is not by any means restrictive since 
the SDE (1) can be transformed into one of unit diffusion coefficient for the 
process Y = {Yj; <t<T} defined as 



where z is an arbitrary element of the state space of X. 

2. A general rejection sampling algorithm. Rejection sampling (RS) is 
a widely used simulation technique. It is frequently presented as follows. 
Assume that f,g are probability densities w.r.t. some measure on R'^ and 



Xt+A =Xt + b{Xt)A + a{Xt) ■ Af{0, A) 
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that there exists e > such that < 1. Then the iterative algorithm: 

Rejection Sampling 

1. Sample Y ^ g. 

2. Sample C/~Unif (0,1). 
3i. If U <e^{Y) return Y. 
32- Else go to 1. 

returns an observation distributed according to /. 

This traditional rejection sampling algorithm seems to imply a fixed order 
in the acquisition of the required random elements: specifically, the proposed 
variate Y precedes the decision variate U. However, such a prescribed order 
for the simulation steps is not necessary. The algorithm to be presented 
in Section 3 is more easily understood when imagining that the proposed 
variate actually succeeds the decision variate. Moreover, this reordering of 
the random steps can be carried out without adding any complexity to 
the algorithm. See [9] for where a similar reordering of the random 

inputs is used to perform an otherwise impossible MCMC algorithm for the 
Dirichlet mixture model. 

Additionally, in some cases there are ways of constructing a condition 
for the acceptance or the rejection of the current proposed element Y from 
minimal information about it. This will be essential in the diffusion context 
where it will never be possible to store a complete, continuous path. A similar 
idea appears in [4] on a perfect simulation algorithm of point processes. In 
that case a probability related with a complicated surface is expressed as the 
probability of an appropriately constructed event whose truth or otherwise 
is easy to verify. 

We now present a formal definition of the RS algorithm in a way that 
incorporates the observations just described. Let {S,S) be a sufficiently reg- 
ular measurable space and i', /x probability measures on it such that is 
absolutely continuous w.r.t. u. Assume that there exists e > such that 
/ := < 1 v-a.s. and that it is easy to sample from u. The following 
proposition can be used to return draws from fi. 

Proposition 1 (Rejection sampling). Let {Yn, In)n>i be a sequence of 
i.i.d. random elements taking values in S x {0,1} such that Yi ^ u and 
P[/i = l|yi = y] = f{y) for all y G S. Define r = min{i > 1 : /j = 1}. Then 
P[Yr£dy] = f,{dy). 

For the proof see the Appendix. 

This presentation of the RS scheme does not assume any order for the 
simulation of Y and I and, besides the certain conditional property given in 
the proposition, does not restrict in any other way the construction of the 
binary indicator I. 



4 



A. BESKOS AND G. O. ROBERTS 



3. The Exact Algorithm. Consider the stochastic process X = {Xt;0 < 
t < T} determined as the unique solution of the SDE 

(2) dXt = a{Xt) dt + dBt, 0<t<T,Xo = 0. 

The drift function a : R i— > R is presumed to satisfy the regularity properties 
that guarantee the existence of a global, weakly unique solution for (2). In 
particular, it suffices that a is locally Lipschitz; that is, for each M > there 
exists Km > such that 

\a{y) - a{x)\ < KmIv - x\; \x\ < M, \y\ < M 

with a linear growth bound; there exists K > such that 

Hx)]"^ <K^{l + x^), xeK. 

See Chapter 4 of [5] for a detailed presentation of weaker conditions. 

Since we are going to carry out rejection sampling it is convenient to think 
of the stochastic processes involved as measures induced on the space C of 
continuous functions from [0, T] to R. We denote by a; a typical element of 
this space. Consider the coordinate functions Bt{uj) =io{t),tG [0, T], and the 
cr-field C = a{{Bt;0 <t< T}). To avoid confusion between the coordinate 
functions and the Brownian motion process, we use the generic notation 
BM = {BMt;0 < f < T} for a Brownian motion (BM) started at 0. Let W 
be the Wiener measure on (C,C) so that B = {Bt; <t < T} is a BM. 

3.1. Rejection sampling for SDEs. We explain how a rejection sampling 
algorithm can be set up for the case of the SDE given in (2). The analysis be- 
gins with the Girsanov transformation of measures. Let Q be the probability 
measure induced on (C,C) hy X = {Xt,^ <t<T}. 

Proposition 2 (Girsanov transformation). Assume that the drift coef- 
ficient a satisfies Novikov's condition: 

Ew exp'^ i f a^{Bt )dt)- < oo. 

It is then true that 



expji J a^{Bt)dt^ 



(3) ^(^)=expy\iBt)dBt - \ j\\Bt)dt^ =: G{B). 

Proof. See, for instance, [8], Chapter 8. □ 

Our goal is to implement a rejection sampling algorithm using (3) to 
construct an accept-reject mechanism. The difficulty is that exact evaluation 
of G{B) is impossible. We can simplify G{B) using Ito's formula to remove 
the Ito integral term and we shall see that this allows us to carry out the 
rejection scheme indirectly. We now need our first assumption. 
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Under Condition 1, G{B) admits the following simplification: let A{u) = 
a{y) dy, « G R. Ito's formula then gives that 

r a{Bt)dBt = A{BT)-A{Bo)-l a'{Bt)dt. 
Jo Jo 

We can now write G{B) as 

G{B) = exp!^A{BT) - A(So) - h + «'(^*)) ^*}- 

Rejection sampling using Brownian candidates is only conceivably possible if 
G{B) is almost surely bounded and this is likely to require A to be bounded. 
To remove this requirement we introduce a third probability measure which 
will be used to construct the candidates for the rejection sampling scheme. 

Consider the biased Brownian motion BM = {BMt;0 <t< T} heuristi- 
cally defined as {BM\BMt = p) with p distributed according to some density 
function /i:Ri-^ [0, oo) w.r.t. the Lebesgue measure. 

Proposition 3 (Biased Brownian motion). Let Z he the probability 
measure induced by BM on (C, C). // the support of h is the real line, then 
Z is equivalent to W and 

dZ , , h(BT) 



' {l/V2^)eM-B^/{'^T))' 



For the proof see the Appendix. 
It is now trivial that 



dQ. , dQdM. . 
dZ^ ' dM dZ^ ' 



where oc implies that we omitted some factors not depending on uj. 

Condition 2. exp{ A{u) — u"^ /2T} du =: c< oo. 

Under Condition 2 and after choosing h{u) = exp{A(ti) — u'^ /2T}/c, 

r-r/i 



(4) iH«-p{-/„ ir'(B,)^yiB. 



dt 



Assume now that the functional involved in the above integral is bounded: 
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Condition 3. There exist constants fei, /c2 S R such that ki < ^q^{u) + 
^a'(ti) < k2 for any u G R. 

We can then write (4) as 

(5) ^(^)ocexp|- 

for a function 4>> defined as 4>{u) = ^a^(u) + ^a'{u) — ki, n G R. This 
creates the possibihty of performing rejection samphng with candidates from 
Z in order to sample from Q, which is the objective. We can choose the length 
r > of the time interval under consideration so that 

(6) < 0(m) < T'^ for any u G R. 

Just consider any T < 1/(^2 — ki) = 1/R where R is the identifiable range 
of (a^ + a')/2 in the sense that the lower (upper) bound we can actually 
obtain for (a^ +a')/2 can be less (greater) than its maximum (minimum) 
value. From this point on it is assumed that T has been fixed so that (6) is 
true. 

We remark that Condition 3 is trivially implied by 

Condition 3'. There exist constants k[,k2 G R such that k[ < a{u), 
a'iu) < k'2 for any u G R. 

Furthermore, Condition 3' implies that Condition 2 holds and moreover 
the constants k'l and ^2 point the way to implementing simple rejection 
sampling algorithms for simulating the endpoint of the biased Brownian 
motion. 

3.2. Constructing the Exact Algorithm. We set H{u!) = Jq (j){Bt) dt. We 
assume from now on that we have access to paths a; Z of the biased 
Brownian motion BM . Note that we can realize such a path at any finite 
collection of time instances by drawing first its ending point ujt ~ h and then 
the rest of the skeleton according to the dynamics of a Brownian bridge. 
Drawing from the univariate distribution h cannot be a big problem; [1] 
gives many algorithms for drawing from densities on R that could be used 
as envelopes for a rejection sampling scheme on h. 

The preliminaries of Section 3.1 together with the general rejection sam- 
pling protocol of Proposition 1 ensure that the following algorithm (were it 
implementable in practice) would output realizations of the diffusion X that 
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solves the SDE (2): 
Impossible Algorithm 

1. Sample a complete, continuous path of BM, w ~Z. 

2. Compute i?(a;). 

3. Produce a binary indicator / such that P[/ = 1\uj] = exp{— iJ(a;)}. 

4. If / = GO to 1. 

5. Output lo. 

The indicator / is easily constructed with the use of some U ~ Unif(0, 1). 
In practice, we can only simulate the path u; at any given finite collection of 
instances <ti,t2, ■■■ ,tn so evaluating the integral H(uj) is impossible. 
However, given that Conditions 1-3 are satisfied, we can produce an algo- 
rithm which manages to circumvent steps 1 and 2 and still carry out steps 
3 and 4 exactly given only a finite but random skeleton of instances of to . 

The idea builds on the simple observation that for a bounded function 
< (p{u) < events of probability (j){u)du can be constructed simply 
by drawing a random point (y, M^) ~ Unif[(0,T) x (0,T~^)]. Then the event 
^ will have the required probability. 

To extend this idea to an event of probability exp{—H) we exploit a 
Taylor series expansion construction which gives us an event of probability 
exp{—H) as an event which depends upon a countable sequence of events of 
probability H. Furthermore, it turns out to be possible to express this event 
both as the countable union of a sequence of increasing events and as the 
countable intersection of another sequence of decreasing events. Crucially, 
for each event in either of the two sequences its truth or otherwise can be 
confirmed by a finite skeleton of a path a; ~ Z and, consequently, the truth 
or otherwise of the event of probability exp{— ff(cj)} can also be determined 
after finite computations. 

All the above ideas are presented in a rigorous way in Theorem 1 that 
follows. The construction to be described is similar in spirit, though in a 
different context, with Von Neumann's comparison method for the simula- 
tion of exponential random variables; see [2] for a detailed review. We have 
denoted by (fi,^, Prob) the underlying probability space that generates all 
the random elements involved in the theorem. 



Theorem 1. Let lo ^Ij he a path of the biased Brownian motion BM on 
[0,r]. Let T = {Vi,Wi)i>i be a sequence of i.i.d. points uniformly distributed 
on (0,T) X (0, 1/T). Consider also some ?7 ~ Unif(0, 1). Assume that lu, t 
and U are independent. Define the following events: 

To = n, r„ = {(l>{Bv, {uj))>Wi,..., HBvA^)) >Wn,u<-. 
(7) ^ ^' 

n = 1,2,.... 
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Consider the sequence of events {En)n>i defined as 

E2n+i = {To - Ti) + (Fa - Ta) + • • • + (Fsn - Tsn+i), 

n = 0,l,..., 

(8) 

E2n+2 — — Ti) + — Fs) H h (r2n — r2n+l) + ^2n+2, 

71 = 0,1,..., 

where (+) implies union of disjoint sets and D — F = D f] F'^ for any sets 
FCD. Then: 

(i) iE2n+i)n>o, {E2n+2)n>o o-i"^ Sequences of increasing and decreasing 
events, respectively, with E2H+1 Q E2X+2 for any k, A G {0,1,...} and 

Prob[nS" E2n+2 - E2n+l] = 0. 

(ii) Let E = E2n+i ■ If I is a binary indicator such that 1 = 1 when E 
occurs and otherwise, then 

Prob[/ = l|a;] = expj- jj^"^ dtj. 
See the proof in the Appendix. 

Recall that Bt is the coordinate mapping Btiuj) = uj{t), t G [0,T]. We will 
from now on write u:{t) instead of Bt{uj). In practice, we can simply iden- 
tify E2n+2 - E2n+i with the null event so that E = [j'^ ^2n+i = 
ricf -^2n+2- Figure 1 illustrates the decomposition of 17 over the sets Ei, i>l. 

Theorem 1 does not impose any restriction on the order of the realization 
of the random elements when rejection sampling will take place. To trans- 
form the theorem into a feasible rejection sampling algorithm it is necessary 




Fig. 1. An illustration of the sets Ei,i > 1. When some even-numbered Ei does not occur 
or an odd-numbered Ei occurs, E does not occur or occurs, respectively. 
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that we draw U and then generate the path lu Z and the sequence r in 
parallel. We follow an iterative process that involves drawing {Vi,Wi) uni- 
formly from (0,T) x (0, 1/T) and then simulating ioiVi) conditionally on the 
already obtained a;(Vi), a;(V2)> • • • ,^iVi-i), for i>l. Recall that this recur- 
sive construction of the path of BM is straightforward as long as someone 
begins by simulating its ending point uj{T) ~ h (the choice for h is given after 
Condition 2 of Section 3.1). Then, the rest of the path is a Brownian bridge; 
given that the locations {cj(Vi), . . . , cj(r)} have been constructed 
the path can be realized at the instance Vi just by drawing 

nn Kr( ni \ ^ ^(^+) - ^(^-) ... t/ ^ 
.:iyi)^M\^{V^) + ^^—^^{v. - j , 

where we have defined VL = max{0, Vj , j = 1, . . . , i — 1 : Vj < Vt} and V- = 
min{T, Vj , j = 1, . . . , i — 1 : Vj > Vi}. See, for instance, page 360 of [3] for the 
representation of a Brownian bridge as a transformation of an unconditional 
Brownian motion which implies the above formula. 

From the definition of (£'i)i>i it is clear that after j iterations are car- 
ried out we have the necessary information to decide if any of the events 
Ei,E2, ■ ■ ■ ,Ej occurred or not. We perform iterations until the first time 
that an odd-numbered occurs or an even-numbered Ei does not occur. 
It is then clear from Theorem 1 that the former case gives 1 = 1 and the 
latter 7 = 0. Further realization of the random elements will not change the 
decision about /. If needed, we can continue simulating an accepted path 
of BM so that the resulted path of X is constructed at any time instances 
requested. 

The recursive definition of {Ei)i>i allows for a simple way of carrying 
out the successive steps of the iterative method described above. Assume 
for instance that 2n iterations have taken place without a decision about 
/, that is, E2n happened and E2n-i did not happen. From (7) and (8) it is 
clear that 

E2n = E2n-1 + ^2n and 
E2n+1 = E2n — ^2n+l = -E'2n-1 + (r2n — r2n+l)- 

The first equation indicates that T2n occurred and the second that we only 
need to check if the subevent of T2n 

T2n - r2n+l = T2n H | V'2n+1 )) < M^2n+1 OI U > 

occurred [i.e., (/)(u;(V2n+i)) < W^2n+i or U > (2n+i)! ] °^ *° reach to a 
similar conclusion about £'2n+i- The same convenient interpretation appears 
for the case when an even-numbered iteration, say the 2nth, is carried out. 
Given that / is not determined before that step we can conclude that E2n 
did not take place if C/ > l/(2?i)! or 4>{io{V2n)) < W2n- 
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We can now present the pseudocode that implements Theorem 1 to carry 
out rejection samphng: 

Exact Algorithm 

1. Initiate a path of BM: Set lj(0) = and draw uj{T) h. 

2. Draw U ~ Unif (0, 1) . Set i = 0. 

3. Draw Unif [(0,T) X (0,1/r)]. SETi = i + l. 

4. Construct mCV) given the currently simulated instances 

OF UJ. 

5i. If (l){uj{V)) <W orU > then 
If i IS even set / = and go to 1. 

If i IS ODD set 1=1 AND GO TO 6. 

52- Else go to 3. 

6. Output the currently simulated instances of uj. 

As already mentioned, we can interpose a step that constructs the proposed 
path at any time instances we require. It is now clear that the algorithm 
returns exact skeletons of the target process X at any given finite collection 
of time instances. 

3.3. A factorization for Q. From a probabilistic point of view, the Exact 
Algorithm manages to decompose the target measure Q in terms of a product 
of Brownian bridges after appropriately extending the underlying probability 
space. 

Analytically, let S be the random skeleton produced by the Exact Algo- 
rithm [we do not include in S the starting point (0,0)]. If Sk is the space 
of the possible configurations of k points {xi,X2, • • • , x^} C [0, T] x R, /c > 1, 
then S takes values in Sk- We denote by s = {{ui,yi), (w2, 2/2), • • • , {uk^Uk)} 
a typical element of this state space. We avoid details about the cr-algebra 
construction on [j^Sn, and simply denote by £5 the distribution of S. Let 
BB(s, x; t,y), for < s <t, x,y £ R, be the probability measure correspond- 
ing to a Brownian bridge starting at the time instance s from x and finishing 
at the time instance t at y. In terms of probability measures, the Exact Al- 
gorithm manages to factorize Q in the following way: 

A: 

(9) q = I^MM{ui-i,yi.i;ui,yi) <g) Csids), 

1=1 

where {uQ,yQ) = (0,0). Critically, the rejection sampling construction of the 
Exact Algorithm allows for the simulation of Cs, so drawing from Q is then 
straightforward. Figure 2 gives a graphical illustration of the factorization 
(9). Because of this decomposition of Q it is possible to identify characteris- 
tics of the target process X after considering the properties of the Brownian 
bridges that fill in its skeletons. Thus, it is possible to simulate exactly hit- 
ting times, extremes (we present these applications analytically in Section 



EXACT SIMULATION OF DIFFUSIONS 



11 



6) and any other random elements for which there are exphcit results for 
the case of Brownian paths. 

Another simple example when we can exploit (9) is at the Monte Carlo 
evaluation of the expected value of functionals of Xt for some time instance 
t € [0,T]. Assume that {Xj,S*}i<j<n are n draws for the skeleton S and 
Xt produced by the Exact Algorithm. An estimator of E[/(Xj)], for some 
function /, can be the mean of the f{Xl), 1 <i <n. Using the simple con- 
ditional expectation property Var[/(Xt)] > Var[[E[/(Xj)|S]] we can produce 
an unbiased estimator of smaller variance after considering the mean of the 
E[/(A't)|S*], 1 <i <n. Conditionally on the skeleton, Xt follows a normal 
distribution, so for reasonable functions / finding the E[/(A't)|S] will be 
straightforward. 



4. Efficiency of the Exact Algorithm. Two aspects of the Exact Algo- 
rithm need to be examined. The first involves the probability of accepting a 
proposed path of BM as a path of X . The second has to do with the number 
of points in (0,T) x (0, 1/T) required to reach a decision about a proposed 
path. 

We can rewrite (5) as 

e(T)^M=exp{-^%W,d*} 

for some appropriate e(T). Note that Q and Z are both probability measures 
so e(T) equals precisely the probability that a path u; ~ Z is accepted as a 
path from Q (see the proof of Proposition 1). Equivalently, e(r) = Prob[/ = 
1] for the binary indicator / defined in Theorem 1. 



X 



X 



X Skeleton S 
< — * Brownian Bridges 



X 

(«2'W2) 



Fig. 2. The factorization ofQ: drawing from Q is achieved after simulating the skeleton 
S = (112,2/2), • • • , (wfc, j/fe)} (in the case of the figure, k = 4:) and then filling in the 

rest of the path with independent Brownian bridges. 
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Proposition 4. For any appropriate T, the probability £{T) of the event 
{I = 1} for I defined in Theorem 1 is at least e~^ and is decreasing in T 
with lim'rio^(^) = 1- 

For the proof see the Appendix. 

As expected, the probabihty of accepting the proposed paths of BM in- 
creases when we decrease the length of the time period under consideration. 

It is important that we give a result for the number of points uniformly 
drawn in (0, T) x (0, 1/T) required to accept or reject a proposed path. Recall 
that Prob) is the underlying probability space that generates all the 

random elements involved in Theorem 1. 

Proposition 5. Define the random variable N {1,2, . . .} as 

iV(?i;) = = - ■''^^^n}, ifw^E", 

\ min{n = 1, 3, G £"„}, ifw^E. 

Then E[iV] < e. 

Proof. Just note that Prob[A^ > 7i\ < Prob[C/ < (^^^\y ], so it is straight- 
forward that 

oo -, 

E[iV] < y 7 TT = e. 

^ n-1! □ 

The random variable N counts the number of events we have to con- 
sider before deciding if w belongs to E or E'^ (equivalently, if the corre- 
sponding realizations of the random elements involved in Theorem 1 make 
E happen or not). It is a perhaps surprising result that for any eligible T 
and any drift coefficient a we are expecting on average less than three of 
the points drawn uniformly from (0,T) x (0, 1/T) before we decide for the 
acceptance or the rejection of a proposed path. 

Note that the Exact Algorithm can take advantage of the Markov property 
of the process X and produce skeletons of any requested length / > after 
merging skeletons of lengths acceptable by the algorithm. The choice of the 
length of the merged skeletons and, subsequently, the efficiency of the Exact 
Algorithm depend on the identifiable range of the functional (c? + a')/2 of 
the drift function. Recall that we can identify analytically which satisfies 

sup (a^ + a'){x) /2 - inf (a^ + a'){x) /2 < R. 

The following proposition gives a result for the case when the Exact Algo- 
rithm merges skeletons of the maximum eligible length T =l/Rto obtain a 
skeleton of length I. \u\ is the minimum integer not smaller than u S R. 
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Proposition 6. If R is the identifiable range of the functional (a^ + 
a')/2 of the drift a and Ni is the total number of the uniformly drawn points 
needed for the Exact Algorithm to return a skeleton of length / > 0, then 

E[Ni] <\l-R\x e^. 
For the proof see the Appendix. 

In total, the Exact Algorithm requires elementary computational skills. 
Except for the appealing characteristic of being exact it seems to compete 
with conventional approximation techniques even in terms of time efficiency. 
The example that follows favors this assertion. 

5. Applying the Exact Algorithm. We apply the Exact Algorithm to the 
SDE: 

(10) dXt = sm{Xt) dt + dBt, < i < T, Xq = 0. 

Conventional methods can only approximate sample paths for the solution X 
of (10) after resorting to one of the suggested time discretization techniques; 
the Exact Algorithm returns exact skeletons of X. 

The drift coefficient a = sin satisfies Conditions 1 and 2. It is also easy 
to check that — 1/2 < isin^(n) + icos(n) < 5/8 for any n G R, so Con- 
dition 3 is satisfied for ki = —1/2 and /c2 = 5/8. In the present context 
(j){u) = ^sin^(u) + icos(n) + ^. We can now choose T = l/{k2 — ki) = 8/9; 
for this ending time instance, < (p{u) < for any real u. The BM 
process is defined w.r.t. a BM via the rule BM = {BM\BMt = p) with 
h(x exp{— cos(n) — n^/2T}. We can draw efficiently from this univariate 
distribution using rejection sampling with Gaussian proposals. 




Fig. 3. A case when the Exact Algorithm accepts the proposed path uj ^1. 
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i« I *■ 



0.* 




Fig. 4. A cose ui/ien the Exact Algorithm rejects a proposed path uj ' 



Everything is now set up for applying the Exact Algorithm. We run the 
algorithm until we generate 5000 exact skeletons of X. We had to pro- 
pose 12,320 paths of BM to get the exact paths so we can estimate that 
Prob[/ = 1] for the indicator of Theorem 1 is close to 0.41. In 58% of the 
proposed paths the decision about accepting or rejecting a path was taken 
after simulating 1 and 2 points respectively uniformly from (0,T) x (0,T~^). 
The maximum number of points needed for the acceptance or the rejection 
of a path was 7 and 6, respectively. 

Figure 3 shows on the left an exact skeleton of X and on the right the 
same skeleton after considering the transformation x i— > </>(a;) for each of its 
joints. In the case of the (/)-path the square encloses the area (0, T) x (0, T^^) 
and the black spots show the location of three points uniformly drawn from 
this rectangle. The numbers next to each circle show the order with which 
they were obtained. At the top right corner we have written the draw U ^ 
Unif(0, 1) needed by the Exact Algorithm. The third point exceeded the 
graph of (j) so the algorithm decided that the event £"3 of Theorem 1 occurred 
after realizing the proposed path lo only at three time instances. The rest of 
the path of X involves the time instances O.Oli for all z = 1, 2, . . . such that 
O.Oli < T and is produced after filling in the proposed path. 

Figure 4 shows similar graphs for the case when a proposed path is re- 
jected. The graph that involves 4) indicates that the proposed path was 
rejected because the event £'4 did not occur. 

Figure 5 shows the estimated density for the distribution of Xi from 
samples of size 1,000,000 obtained after using the Exact Algorithm and the 
Euler approximation method for different time discretizations. For the case 
of the Exact Algorithm, to draw from Xi we had to merge skeletons on 
the time intervals [0,8/9] and [8/9, 1]. Table 1 presents the times in seconds 
needed to get these samples and the p-values of the Kolmogorov-Smirnov 
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test that compares the approximate draws of the Euler method with the 
exact ones of our algorithm. All programs were written in C-language and 
executed on a Athlon PC, running at 1500 MHz. It is clear, at least for the 
case of the SDE given in (10), that the Exact Algorithm is superior to the 
Euler approximation even in terms of computational time. 

Note also that, as implied in Proposition 6, the time needed to draw from 
Xi, I > 0, increases linearly in I for the Exact Algorithm. In contrast, the 
Euler method needs thinner approximation to produce reliable results as / 
increases because of the accumulating errors. 

6. Exact simulation of extremes and hitting times. We take advantage of 
the properties of the Exact Algorithm to simulate exactly the maximum and 
the hitting time of a horizontal line boundary for one-dimensional diffusions 
determined by SDEs of the type (2) under the Conditions 1-3 given in 
Section 3.1. 

As explained in Section 3.3, the Exact Algorithm reduces the problem of 
detecting characteristics for a path of the process X to the much more 
straightforward task of carrying out the certain detection process for its 
Brownian subpaths. To simplify the presentation of the algorithms that fol- 
low we denote by Si{x;ti,t2, ■ ■ ■ ,tn} a skeleton of the target process X at 
the time instances < ti < t2 < ■ ■ ■ < tn = I starting at x. 




-A -2 2 4 6 



Fig. 5. The estimated density for the distribution Xi from samples of size 1,000,000 
generated by the Euler approximation (four cases for time increments , , 2^^, 
2^^ ) and the Exact Algorithm. 
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Table 1 

Times in seconds needed for the Exact Algorithm and the Euler 
approximation to produce 1,000,000 draws from the distribution of Xi 



Exact 


Euler 

(for inc. A = 2"" , n = 2, 3, . . 


.,8) 


KS test 




5 sec 









11 sec 









15 sec 







35 sec 


24 sec 









45 sec 




0.002 




87 sec 




0.261 




174 sec 




0.412 



The last column shows the p-values for the Kolmogorov-Smirriov test with 
null hypothesis that the exact and the corresponding approximate draws 
come from the same distribution. 



For the case of the maximum value of X over the time interval [0,/], 
/ > 0, we apply the Exact Algorithm, if necessary after dividing [0, /] in 
smaller pieces of length at most the maximum length T permitted by the 
Exact Algorithm and dealing with each piece separately, and locate the 
maximum value of the accepted proposed paths of BM after drawing the 
maximum of all the Brownian bridges (BBs) that intervene between the 
successive unveiled instances of these paths. Drawing the maximum of a BB 
is straightforward. 

Let BM^ = {BM^'jO < s < t} be a Brownian motion over the interval [0, t] 
started at y. If Mt = sup{ i?M|; s G [0,t]}, then it is a known result that 

P[Mt e dblBMf = a] oc (26 - y - a) exp|- ^^^~|^~ "^^ | db, 

b > max{y, a}. 

For a proof, see, for instance, [3], page 95. This is just a linear transformation 
of a Rayleigh distribution and it is easy to verify that 

[Mt\BM^ = a] = l(V2tE(l) + (a-y)2 + a + y), 

where E(l) denotes an exponential random variable with unit mean. This 
formula generates the maximum of a BB of length t starting at y and fin- 
ishing at a (for arbitrary t, y and a) and will be applied for all the BBs that 
fill in the exact skeleton of X. Thus, the complete algorithm for drawing the 
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maximum of X = {Xt;0 < t < 1} is as follows: 

Exact Simulation of the Maximum 

1. Call on the Exact Algorithm and draw a skeleton 
Si{0;ti,t2, ■ . ■ ,tn}- 

2. Simulate the maxima M^,M^,...,M" of the BBs that fill in Si. 

3. Output sup{M*; i = 1, 2, . . . , n}. 

A similar procedure is carried out for the simulation of first passage times. 
We now have to check if each of the BBs that fill in the skeleton of X 
produced by the Exact Algorithm hit some arbitrary boundary 7 or not and 
for the first bridge that hits 7 to find the precise time instance when that 
occurs. Assume for simplicity that 7 > is greater than the starting point 
of X. 

The hitting time of a horizontal boundary for a BB is closely related 
with the first passage of an unconditional Brownian motion over a linear 
boundary. It is a known result that if BM (6) = {BMs{5);0 < s < t} is a BB 
of length t started at and finishing at 5 and BM = {BMg; s > 0} is an 
unconditional Brownian motion started at 0, then 

(11) BMs{6)lU+*-^BM,/^t_,), 0<s<t. 

Let T^{6) = inf{s G [0, t] : BMs{5) > 7} and r^,^ = inf{s > : BMs >ri + Cs} 
with the convention that inf = 00. Then for any s G [0,t] it is true that 

P[Ty{6) >s]= F[BMu{6) < 7, for ah < n < s] 

u 



u ^ t 

t 



BM^* < ^ + ^ 



BM^I(^i_u) < 7) for all < u < s 
6 



-u 



for all < u* < 



t 



> 



for 1] = "f/Vi, C = (7 ~ ^)/Vt- It is now clear that for these values of rj, 



(12) T^(5)^5(r,,c) 



for 5(n) = |*^/(^ + ^)' 



<u < 00, 
u = 00. 

The density of r^^^ is given by the Bachelier-Levy formula: 



(13) 



p'^'Hn) 



exp{-(Cn + 7?)''/2n}, u > 0. 



A proof can be found in [6], Chapter 1. If we denote by IG{^,X), /x > 0, 
A > 0, the inverse Gaussian distribution with density: 



IG{fi,X,u) 



A 



2ttv? 



exp 



u>0, 
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then, when r]( < the density (13) is easily identified as an IG{—r]/C,,rf') 
and integrates to 1, that is, the Brownian path hits the hnear boundary 
rj -\- Qt with probabiUty 1, as expected. When 77C > 0, (13) can be written 
as exp(—2(^r])IG{r]/(^,r]'^,u), which means that the Brownian path hits the 
Hnear boundary w.p. exp(— 2(^r/) and when it does the distribution of the 
hitting time is just the IG{ri/C,^rf'). Draws from the inverse Gaussian distri- 
bution can be generated in a very efficient way; see the algorithm described 
in [1], Chapter IV. Thus, simulating r^^^ and, via (12), t^{5) for any values 
of the involved parameters is straightforward. 

We now have all that is needed to carry out the algorithm for the exact 
simulation of the hitting time of a horizontal boundary for the solution X 
of (2). Since it is not possible to know a priori the length of the path of X 
we need to construct before we locate the time instance when X hits the 
boundary, we merge as many exact skeletons of X of eligible length T as 
necessary before some intervening BB hits the boundary. 

Exact Simulation of Hitting Time 

1. Set i = l, y = 0. 

2. Call on the Exact Algorithm and generate 

31. If the intervening BBs hit 7, save the first time t when 

THAT OCCURS. 

32. Else, set y = Xt, i = i + l and go to 2. 
4. Return (i - 1)T + t. 

Assume that N'^ is the total number of the uniformly drawn points needed 
for finding = inf{t > : = 7} and that Ni is the number of these points 
needed for accepting a proposed path from (i — 1)T to iT, i > 1. Let Gi be 
the information for the obtained exact path of X until the time iT, i>0. 
Proposition 6 is true for any given starting point of the target process, so 
E[iVj|^j_i] < [T • R\ for all i > 1, for R as defined in Proposition 6. We set 
T := \t^/T~\ ; r is the total number of the exact skeletons we have to merge 
before we find r^. Clearly, 



00 00 
E[N^] = E[Ni ■ I{t > i}] = E E[E[iVi • I{r > ijlGi-i]] 

i=l 1=1 

00 00 

= J2 HHr > i} ■ E[iV,|g,„i]] <e^\T-R\Y,nr>i] 

i=l 1=1 



<e^rr.fll(5M + i), 

where we have used the fact that {r > i} is ^j„i-measurable. So the expected 
time for the termination of the algorithm, in terms of the uniformly drawn 
points needed, is finite when E[r^] is finite. 
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6.1. An application. We applied the above Exact Algorithms to the solu- 
tion X of the SDE (10) considered in Section 5. We compare the exact draws 
of our algorithms with the approximate ones of the simple Euler scheme that 
considers the continuous time process Y = {Yt}t>Q, Yq = 0, defined on the 
instances {i/i}i>i, for some chosen increment h> 0, via the recursion 

Yih = + sin(y(j„i);j)/i + Zi^h, 

where Zi^^i ^ ^ are i.i.d. draws from the normal distribution with mean 
and variance h. The paths of Y become continuous after considering the 
linear interpolations between the successive instances of the grid {i/i}i>i. 

On the left of Figure 6, we show a ^^-plot comparing two samples each of 
size 50,000 from the distribution of M2' = sup{Xt;0 <t <2} generated by 
the Exact Algorithm and the Euler approximation. For the Euler scheme we 
used increments h = 2~^. On the right of the same figure, we show the times 
in seconds needed to get 50,000 draws from M2' for the Exact Algorithm 
and the Euler scheme for different increments h. We have also included 
the results of a Kolmogorov-Smirnov test that compares the approximate 
samples of the Euler scheme with the exact sample as an indication of how 
small h has to be for the Euler scheme to give correct results. 

In Figure 7 we present similar results for the case of the simulation of the 
first passage time T2 = inf{i >0:Xt = 2}. Actually we tune the algorithm 
to obtain draws from the min{T2, 10} since it can be shown that E[t2] = 00 
so the expected number of the uniformly drawn points for the completion 
of the algorithm for drawing from T2 is infinite. 

It is remarkable that in both cases the Exact Algorithm appears much 
more efficient than the Euler approximation. 



Si 




Exact 







Euler 




Exact 


h 




Times 


KS 


T 


Times 


2" 


-4 


1.2 s 









2" 


-5 


2.4 s 





0.25 


2.7 s 


2" 


-6 


5.0 s 









2-^ 


-7 


7.9 s 





0.50 


3.3 s 


2- 


'8 


17.4 s 









2" 


^9 


32.7 s 





1.00 


3.1 s 


2-^ 


-10 


72.1 s 


0.013 






2- 


-11 


125.2 s 


0.017 






2' 


-12 


245.9 s 


0.006 







Fig. 



solution of (10). 



Results from the simulation of the maximum = sup{Xt; <t<2} of the 
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Euler 



Exact 



Exact 



h 




Times 


KS 


T 


Times 




-5 


5.2 s 









2- 


-6 


11.4 s 





0.25 


9.1 s 


2- 


-7 


22.3 s 


0.00001 






2- 


-8 


43.8 s 


0.00308 


0.50 


5.9 s 


2- 


-9 


84.1 s 


0.01404 






2- 


^0 


170.6 s 


0.51623 


8/9 


5.1 s 


2- 


-11 


345.4 s 


0.92642 







Fig. 7. Similar results to those of Figure 6 for the case of the simulation o/ mm{r2, 10} 
for T2=mi{t>0:Xt = 2}. 



7. Extensions and conclusions. In this paper we have introduced a simple 
but computationally effective way of simulating exactly from a family of 
diffusion processes. The algorithm outputs a skeleton which can be readily 
"filled in" as and when necessary using simple Gaussian random variables, 
and crucially independently of the diffusion we are attempting to simulate 
from. 

We have not carried out an extensive simulations study. However, in the 
examples considered, our method performs very favorably in comparison to 
the obvious numerical approximation alternative using the Euler scheme. In 
Section 4 we also give results which show that the method can be robust to 
the length of the time-series, and computing time is at worst linear in the 
degree of nonlinearity (as measured by the range of + a'). 

The convenient form of the output allows the algorithm to be used in a 
number of ways including the construction of reduced variance Monte Carlo 
estimation. In Section 5 we discuss unbiased estimation of boundary hitting 
times and diffusion maxima. We envisage future application in inference for 
stochastic processes and finance. 

The most demanding of the Conditions 1-3 (detailed in Section 3.1) re- 
quired for the algorithm to work, is that the functional + a' of the drift be 
bounded. In most cases + a' is bounded from below but not from above 
so we can still produce (5) for some (j)>0 and hope for a valid rejection 
sampling scheme. Ongoing work is investigating such an approach. 

It is worth remarking that the approach outlined in this paper extends 
routinely to SDKs for non-Markov processes absolutely continuous with re- 
spect to appropriate martingales, and our focus on diffusions has been purely 
for simplicity. Furthermore, it is easy to extend these results some way to- 
ward considering jump diffusions. 
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Proof of Proposition 1 (Rejection sampling). Note that for any 
i = 1, 2, . . . we get 



e. 



P[h = 1] = / m = MYi = vHdy) = f f{y)u{dy) = f e^i{dy) 
Js Js Js 

Trivially, for any F € 5 we have the property 

(14) P[Yr £F]= F[Yr £F,h = 1] +F[Yr £ F\Ii = 0] • P[/i = 0]. 

From the independence among the members of the sequence {Yn,In)n>i it 
is clear that 

(15) P[YreF\Ii = 0]=P[YreF]. 

We can easily find the P [Yj- £ F,Ii = 1] in the following way: 

(16) P[Yr eF,h = l] = jp[h = l\Yi = y]v{dy) = f{y)v{dy) = e^i{F). 

From (14) using (15) and (16) we get 

P[Yr G F] = e/i(F) + (1 - e)P[Y^ G F], 
which yields that P[y^ G F] = ^(F). □ 

Proof of Proposition 3 (Biased Brownian motion). Choose any F G 
C. We will show that Ez[If] = Ew[If/] where we have set 



exp(-S|./(2r)) 



Note that / is cr(i?T)-iiieasurable. From the definition of BM it is clear that 

Z[F\a{BT)] = W[F\a{BT)] =: giBr) W-a.s. 

for some Borel-measurable function (7 : R R. Clearly, 

Ew[If/] = Ew[Ew[lF/k(5T)]] =Ew[/W[F|c7(Sr)]] 

h{u)^j2^ exp(-nV(2r)) [ u \ ( \A 

; o ,, r^.. , q(u)du= / n(u)q(u)du 

R exp(-u2/(2r)) ^2^ 

since w.r.t. W the random variable Bt is distributed according to the normal 
distribution with mean and variance T. It is straightforward that 



Ez[lF] = MEz[lF\(r{BT)]]= / h{u)g{u)du 

JR 

since under Z we know that Bt ~ h. It is clear that E^p^?] = Ew[If/]- D 
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Proof of Theorem 1. (i) From (8) it is straightforward that Ei C 
-E's ^ • • • ^ E2n+i ^ -E'2n+2 for any n G {0,1,...}. From the definition of 

(£'2n+2)n>0 We Can conclude that E2n = S2„+2 + (r2n+l -r2n+2), so -E2n+2 ^ 

E2n C • • • C Clearly, E2K+1 C £'2a+2 for any k, A € {0, 1, . . .}. 

The definition of (r„)„>o implies that for any n > it is true that r„+i C 
r„, Cln^n is a set of zero probability and r2n+2 = E2n+2 — E2n+i- Therefore: 

I P|-E'2n+2 I — I [J-E2n+1 J = P|(i?2n+2 — -E'2n+l) = P|r2n+2- 
\ n / \ n / n n 

Trivially, Un E2n+i ^ fin -E'2n+2 and their difference fin -E'2n+2 - Un -E^2n+i 
has zero probability. 

(ii) Since -E'2n+i T -^^ it is true that 

(17) Prob[/ = l\uo] = Prob[£;|(j] = J[iin^Prob[£;2n+i|w]- 
Recall that E2n+i = Yl^k=a^2k — ^2^+1)- We can now get that 

n 

Vioh[E2n+M = ^(Prob[r2fc|c^] - Prob[r2fc+i|c^]) 

(18) 2n+l 

= ^(-l)^'Prob[rfc|..]. 

k=0 

From the fact that r = {Vn-, VFn)n>i are i.i.d. and r, U and uj are independent 
we conclude that 



(19) Prob[rfc|cj] = Prob 



\{VTO^(|,{BvM)>W^\u:]. 



i=l 



The Prob[(/>(i?v'i (1^)) > Wi\u;] is the probability that given a path cj ~ Z a 
point uniformly selected from the rectangle (0,T) x (0, 1/T) is found below 
the graph {{t,cj){Bt{u)))) U £ [0,T]}. Recall that Bt is just the coordinate 
mapping Bt = Lv{t), t G [0, T]. From (6) and elementary probability theory 
we get 

Froh[cj)iBv,{uj)) > Wk\io] = (l){Bt{uj))dt. 

Jo 

It is clear that Prob[?7 < -^\uj] = so (19) yields 

(20) Prob[rfc|a;] = ^|^^ </)(St(a;)) dtj'. 
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Starting from (17) and using (18), (20) we get 

rT ^ fc 



Prob[I=lM=Jim ^jy^ mio^))dtj 

= exp|-^^0(Bt(w))dt|. □ 



Proof of Proposition 4. Since < < T^^, the probability of ac- 
cepting an arbitrary path ~ Z will be exp{— 4'{Bt) dt} > exp(— 1). 

Recall that T is restricted not to be bigger than some constant Tq{= 
-j^^^^ ) ; see Condition 3 for the definition of ki,k2. To emphasize the in- 
volvement of the time variable in what follows we define p : [0, Tq] x C ^ [0, 1] 
with 

p(T, uj) = expj - (j){Bt{u;)) dtj . 

From result (ii) of Theorem 1 it is clear that 

(21) e{T)= I p{T,Lu)Z{dio) 

Jc 

for all eligible T. Recall that < (p < for any T <Tq. Trivially, p is 
decreasing in T and limxio p{T,uj) = 1, both properties being true Z-a.s. 
The first property yields, after using (21), that e(r) is decreasing in T. 
Also, since limn-toop{l/n,uj) = 1 and p{^,u>) < p(;^, w) for any n > 1 the 
monotone convergence theorem gives 

lim e(l/n) = lim / p(l/n,Lj)'Ij(dLo) = / < lira p(l/n,uj) >Z,(dLu) = 1. 

n—joo n— »oo J(j i \ / In— ♦oo J 

From the monotonicity of e(T) we conclude that \im.xiQe{T) = 1. □ 

Proof of Proposition 6. For the time interval [0,T], for T = l/R, 

the algorithm uses Nj- := Ni + N2 H h Nj points, where Ni is the number 

of points needed to decide about the ith proposed path and J is the number 
of the proposed paths until one is accepted, J = 1, 2, . . . and i = 1, 2, . . . , J. 
We denote by the expected number of points for deciding (in general) 
about a path w ~ Z and by E[A^|A] and E[A^|A'^] the expected number of 
points for deciding about a path w ~ Z given that the path is accepted and 
rejected, respectively. Let e = Prob[j4] be the probability of accepting a path. 
Then 



E[Nt]=E[E[Nt\J]]=E 



■j-i 

Y,nN,\J]+E[Nj\J] 

.i=l 



E[(J - 1) • E[N\A^] + E[N\A]] = {l/e- l)E[iV|A^] + E[iV|^] 

-•E[iVl. 

e 
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From Proposition 4 we have that 1/e < e and from Proposition 5 that E[A^] < 
e, so E[A'^j'] < e^. The same result can be obtained for any \l ■ R\ skeletons 
of length T = 1/R (or less, for the case of the last skeleton when 1/R does 
not divide I) that will merge to construct the complete skeleton of length I. 
□ 
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